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Introduction 

This documents collects the results of my studies and my researches during 
the period between April 2012 and November 2013. During this period I 
wrote my Master’s thesis about CT and SPECT (until October 2012), I 
worked on mixed research group on SPECT imaging, and I made a number 
of seminars (May 2013 - November 2013). 

This work is intended for mathematicians (or physics, or engineers) who 
want to learn some basics of medical imaging systems from a technical (non¬ 
medicine) point of view. It could be a good guide for beginners, even if 
incomplete and, maybe, wrong in some parts. 

Medical imaging systems 


Essential glossary 

• analytical: precise, exact. Analytical methods are based on analytical in¬ 
version formulas. 

• attenuation: physical phenomenon of the deviation of an electromagnetic 
ray from its natural path, due to the contact of a body. 

• CT: X-ray Computed Tomography, a medical imaging procedure that utilizes 
computer-processed X-rays to produce ’slices’ of specific areas of the body. 

• dummy: (also called phantom, we use that name to avoid confusions) 
physical device that simulates a certain organ and used to make experiments 
on the tomography machine, in vitro. 

• electromagnetic waves: is an oscillation which propagates in the space, 
transporting energy but not matter. We distinguish different kind of waves 
by their wavelength, among which: 

— X-rays: wavelength in the range of 0.01 to 10 nanometers. 

— Gamma-rays: wavelength shorter than 1 picometer. 

• ill-conditioned problem: given /( x) = y, y € Y, we have to find x £ X, 
with X and Y vector spaces. This problem is i.c. if an arbitrarily small 
perturbation on y can lead to an arbitrarily large error on x. 

• in silico: experiment made on mathematical simulation of the data (in this 
work we use phantoms) and not using real data. 

• in vitro: experiments made on real data coming from a biological/chemical 
phenomenon reproduced in a test tube (in this work we use a dummies), 
and not using an organism. 



• kernel function: a function K such that K{x, y) = <j>{x — y) with (f> a RBF. 

• nuclear medicine: a branch of medicine involving the application of ra¬ 
dioactive substances in the diagnosis and treatment of disease. 

• over(/under)-determined: a problem with more(/less) equations than 
unknowns. Given /( x) = y, y £ Y, we have to find x £ X, with X and 
Y vector spaces. This problem is o.d. if dim(Y ) > dim(X) and u.d. if 
dim(Y) < dim(X). 

• phantom: an image created to test the efficiency of the algorithms in silico , 
without need of data given by a physical machine. The most used phantom 
in this work is the Shepp-Logan phantom, which is similar to a brain, see 
Figure ??. 

• Randon transform: is the integral transform consisting of the integral of a 
function over straight lines, proposed by the Austrian mathematician Johann 
Radon in 1917. Is simulates mathematically the effect of attenuation of 
some electromagnetic waves passing through different materials at several 
angles, and then it is useful for the CT reconstruction. 

• RBF: Radial Basis Function, functions <j> : C R s —>• R. which depends 
only on the radial component r = ||x||. 

• sampling: times or points at which the data is collected from a continuous 
phenomenon. 

• scattering: erroneous acquisition of some rays made by the tomography 
machine due to attenuation. 

• sinogram: an l x p matrix of data collected by the tomography machine in 
l angular scans. In the case of in silico experiments we call s. the discrete 

Radon transform of the phantom. 

• SPECT: Single-Photon Emission Computed Tomography, is a nuclear medicine 
imaging technique using radioactive tracers that emit Gamma-rays in de¬ 
caying. 
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Chapter 1 

Introduction to Analytic 
Tomography 


Axial Tomography machines provide result of some circular acquisitions. 
The aim of reconstruction is to obtain from projection data an image of 
the slice of the body as seen from above. In the next chapter we discuss 
algorithms of reconstruction from projections and their efficiency in a few 
cases 


• transmission tomography (e.g. CT): an electromagnetic ray passes 
trough the patient and is detected at the exit in order to get a mor¬ 
phological analysis of its interior. 

• emission tomography (e.g. PET, SPECT): a radioactive tracer is 
injected into the patient and detected by the machine in order to make 
an internal functional analysis of the organs. 

• hybrid tomography (e.g. SPECT/CT, SPECT/MRI): two simulta¬ 
neous analysis. 


1.1 Resolution and other technical issues 

1.1.1 Resolution 

The spatial resolution of a gamma camera is a measure of its ability to re¬ 
solve small objects in the field of view. Spatial resolution can also be defined 
as the minimum distance between two points such that they can be pictured 
separately. This means that objects placed at a distance smaller than the 
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Figure 1.1: A scheme of the main differences between transmission and emission 
tomography. 


resolution limit are imaged as a single blurred one. Resolution is therefore 
a crucial parameter that measures the reliability of the gamma camera in a 
specific setting [?]. 

The overall spatial resolution of a gamma camera system ( R s ) depends on 
different factors, both geometrical and physical, and it is usually expressed 
in terms of collimator resolution ( R c ) and intrinsic resolution ( R, t ) . R s is 
typically assessed from the full width at half maximum (FWHM) of the pro¬ 
file of a point-like, or line-like, radiation source. Different methods can be 
used to calculate the FWHM of a point or line spread function. 

The system resolution R s depends on the collimator resolution R c and on 
the intrinsic resolution Ri. Using the convolution mathematical theory[?], 
we obtain R 2 S = R 2 + R 2 , which gives us 

R s = ^R 2 + R 2 (1.1) 

due to the fact that R s will be positive. 

The intrinsic resolution R, L is linked to the properties of the detector and 
electronics. For a given photon energy, Ri could be considered independent 
of the object-to-collimator distance whereas the collimator's resolution de¬ 
pends largely on the geometrical layout and can be expressed as a function 
of a number of parameters: 

• x : distance between the object and the collimator's surface; 

• L: collimator’s hole length; 
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• D: collimator’s hole size; 

• c: crystals thickness; 

• t: thickness of the septa; 

where L, D , c, t are declared by the manufacturer as well as Ri. 

Figure 4.1 schematically shows the geometrical layout of a point source. 



Figure 1.2: Geometrical layout of a point source acquisition by a gamma camera 
equipped with parallel hole collimator. 

To calculate R c , we consider gamma rays coming from a point source P 
(as in Figure 4.1) and particularly rays PA (parallel to the septa) and PB 

(angular limit). Now, since the radiation profile is similar in shape to an 

A A A 

isosceles triangle ( HIB ) and triangles HIB and HFG are similar, then the 

FWHM is about half of the base ( R c ~ AB). 

A A 

Because of the similitude of the triangles PAB and P'A'B and the fact that 
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t <C D(= A'B') we infer that 

P 7 ! 7 : WB 7 = RA:AB 


which is equivalent to: 


L : D = (c + L + x) : R c (1.2) 

Instead of L, L e ff is usually used which is a length that is weighted to take 
the septal penetration into account, 

L eff = L ~ ~ ( L3 ) 

The constant /j, is the linear attenuation coefficient of the material of the 
collimators (usually lead, /j = 2.49mm~ 1 at 140 KeV). Thus from (1.2) and 
(1.3): 

1.1.2 Electromagnetic waves 


An electromagnetic wave is an oscillation which propagates in the space, 
transporting energy but not matter. Its intuitive form is a cosine function 
with amplitude A and frequency u. Its speed in the vacuum is c = 3- 10 8 m/s, 
the speed of light, so if A is the wavelength, then \v = c. The energy carried 
by the wave is E = hu (where h = 6 • 10 -34 Js is the Plank constant), hence 
the energy is inversely proportional to the wavelength: 



(1.5) 


Some examples of electromagnetic waves are the light of the sun, infrared, 
microwaves, radio waves, up to X-rays and gamma rays. 


It is known as X-rays that portion of the electromagnetic spectrum with 
a wavelength in the range of 10 -3 up to 10 1 nm discovered in 1895 by the 
German physicist Wilhelm Conrad Rontgen, who earned for this discover 
the first Nobel price in 1901. X-rays are produced from the transition of 
electrons, and they can penetration through many objects, included the hu¬ 
man body. For this reason X-rays are widely used in medicine to discover 
irregularities into the patient’s body. 
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Similar to X-rays, but shorter in wavelength (then they transport more 
energy for the (1.5)), the gamma rays (sometimes written as 7 -rays) are 
produced by a radioactive decaying 1 . As the X-rays, gamma rays can be 
dangerous and cause burns, cancer, and genetic alteration, so doctors must 
pay attention not to do useless or too hazardous exams, by limiting the 
energy of the waves used and the time of exposition. 


1.1.3 Attenuation 

Every electromagnetic waves passing thought a material mean is subject to 
an attenuation due to diffraction; in every point the ray has a probability 
p to deviate from its natural path 2 . The linear attenuation coefficient // is 
proportional to the probability of deviation while passing in 1 cm of matter; 
p depends on the material, on its density and on the energy carried by the 
wave. The linear attenuation coefficient is proportional to the density p, so 
the mass attenuation coefficient p rn = only depends on the energy and 
the material. 

Table 1.1 shows the mass attenuation coefficients of different materials (wa¬ 
ter, sodium iodide and lead) and energies. 


Photon energy 
(eV) 

Mass attenuation coefficient p m 
(cm 2 /g) 

h 2 o 

Nal 

Pb 

80 

0.179 

3.00 

2.07 

100 

0.168 

1.64 

5.23 

150 

0.149 

0.590 

1.89 

200 

0.136 

0.314 

0.945 

300 

0.118 

0.158 

0.383 

400 

0.106 

0.112 

0.220 

500 

0.097 

0.092 

0.154 

600 

0.089 

0.080 

0.120 

800 

0.079 

0.066 

0.085 

1000 

0.070 

0.058 

0.069 


Table 1.1: Mass attenuation coefficient for several values of the photon energy and 
for different materials. 


1 this is the fundamental difference between X and gamma rays: the first ones come 
from outside the nucleus of the atoms, the second ones from inside. 

2 see [?] Ch 6D, p.80 for a complete and detailed survey about photon attenuation and 
scattering. 
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It is evident from this table that the mass attenuation coefficient is de¬ 
creasing with the growth of the energy or with the decreasing of the density, 
but it is neither inversely proportional to the energy or proportional to the 
density. Another popular way to represent this attenuation is the Hounsfield 
scale 

j-j M l^water 

l-lwater 

which compares the attenuation coefficient of the medium and that of the 
water. 

Computerized tomography, in its classical form, calculates the attenuation 
coefficient of the different areas inside the body from the difference between 
entering and outgoing X-rays, in order to recognize the different organs and 
find possible malfunctions or malformations. 

1.1.4 Scattering 

During the reconstruction we have to face different kind of issues related 
to the attenuation. The most known, called scattering, is the effect of a 
multiple attenuation of the same ray. In fact assume that a ray is diffracted 
in the point (x'i, y \) and deviates from its natural path; assume also that, as 
in Figure 1.3, it deviates again in the point (.T 2 , 2 / 2 ) in the direction of the 
detector (perpendicular in the figure). This means that the machine will see 
many rays “out of place”, but with low energies. 

The most common way to reduce the scattering effect 3 is a Hamming win¬ 
dowing of the incoming data around the known peak of energy of the used 
tracer or ray, as shown in Figure 1.3. 




Figure 1.3: Scattering effect (left) and scattering correction example for 99 m Tc , 
who has a single peak of energy of 140KeV (right). 


3 See [?, Ch IOC p.161] for further explanations. 
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1.1.5 Partial volume effect 

Another typical issue in tomography is the correct reconstruction of the 
higher-valued areas. In fact in such area the reconstruction will result 
smoother and lower-valuer than the original ones. Conversely the lower¬ 
valued areas next to it will have a growth of value near the higher area. 
This is known as partial volume effect, and it is similar to the light bulb 
effect: it is hard to see the exact position of of the bulb since it is on because 
of the light rays next to it. As we will see later the partial volume effect is 
not related to the reconstruction algorithm we use, but only depends from 
resolution. 



Figure 1.4: Partial volume effect example on the characteristic function of an in¬ 
terval (right). 


1.2 Mathematical modeling of problem 

1.2.1 Basic assumptions for modeling 

In order to model correctly the tomography machine we will make some 
assumptions: 

• all rays have no width and we can consider them as straight lines. 

• all rays at angle 9 are parallel. 

• the machine runs l scans at the angles (9 \,..., 9i); the detectors move 
on a circle (typically l = 120,60). 
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• for each scan j € 1the machine gets p linear samples (typically 
p = N, or p = 128 or p = 64). 

• the total time of acquisition is relatively short so is indifferent to run 
an acquisition before or after another one. 

• the reconstructed image has dimension N x N (typically N = 128 or 
N = 64). 


1.2.2 Radon Transform 

The Radon transform, introduced in 1917 by the Austrian mathemati¬ 
cian Johann Karl August Radon, can be physically interpreted as a ray 
passing through different materials, each one of them with its own attenu¬ 
ation coefficient. We know the intensity in which the rays entered in and 
released by the body, so we can find the attenuation coefficients. In fact, 
Beer-Lambert’s law says that if I(x) is the intensity of the ray, A{x) the 
attenuation coefficient at the point x, satisfies 

A I = —A(x)I(x) Ax. 


Equivalently, by integrating their respective differential terms, 


fXl 


A(x)dx 


>x 0 



log 


-fOci) 

I(xo) 


( 1 . 6 ) 


where the second term is known. Equation (1.6) has an equivalent on the 
2-dimensional plane. In fact, given the change of coordinates due to rotation 


x 

V 



— sin 6 \ 
cos 6 J 


t 

s 


we get the formulation of the Radon transform as an integral on a line 


Kf(t, 0) = I f = [ f(x)S(t - x ■ 9) dx 

4R2 

where x = (x,y), 9 = (cos 9, sin 9), or equivalently: 




f(t cos 6 — s sin0, t sin0 + s cos 9) d-s 


f(t.9 + s9 ± )ds (1.7) 


where 9 ^ = (— sin0, cos 9). 
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1.2.3 Backprojection 

Now we ask for every point (x, y ) which is the average of the rays that pass 
through that point. This question is answered by the adjoint operator to 1Z. 
called Backprojection operator 

K*g(x,y) =- [ g(xcos0 + y sin 9,9) dO = 

P I Js 1 

tfu [ g{s,9)5(s - x ■ 9) dsd0 

P I JVLxS 1 

where S 1 = [0,7r] or S 1 = [0 , 27t] 

Warning! 1Z* is not the inverse transform of 1Z. In fact 

n*nf{x) = 7ZT * /. 

PI 

This means that using the backprojection operator with projection data will 
lead to a blurred image. 
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1.2.4 Attenuated transform 

In the case of emission tomography things are more complicated. From 
Beer’s law we obtain that 

I(x\) = I{xo) exp A(x)dx' S j 

Suppose to know the attenuation coefficient, say a(x), we want to obtain 
the radioactivity f(x) by its angular projections. 


We define the attenuated Radon transform (briefly AtRT) 


K a f(t,0)= [ e~ Va ^ e+ ^f{x) 

where D is the Divergent beam transform defined as follows 


Vh(x, 9) 



h(x + t cos 9, y + t sin 9) dt 



h(x + 19) dt 
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Chapter 2 

Reconstruction in X-ray 
computed tomography 


From the introduction given in the previous chapter we know that we have 
to solve two kind of problems: 


1. in the CT case (transmissive) we have to solve the following problem 
RT problem 

Given g projection data 
find / such that 7 Zf = g 


2. in SPECT and PET case (emissive) we can approximate / with the 
solution of the previous problem or solve 
AtRT problem 

Given g projection data and a attenuation map 
find / such that 1Z a f = g 


to estimate the attenuation map we may need a simultaneous CT to¬ 
mography. This is the reason why we may need two simultaneous analysis 
like SPECT/CT or PET/CT. 


This chapter is about the first problem, while the second one will be 
treated in the next chapter. 
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2.1 Analytical methods 

2.1.1 Some theorems 

Proposition 2.1.1. The hackprojection operator is the Hermitian added to 
the transform 77, but not its inverse: in fact 

n*nf(x) = 2 [ rl^-dy = *2 / 

J R 2 \x - y I \x\ 

Proof, see [3] Ch 2.1, p.15 for the proof. □ 

this means that what we get from applying the hackprojection operator 
to projection data is a blurred image. From this property we can derive 
a first numerical method of inversion which simply uses the 2D FFT and 
IFFT algorithms for a fast deconvolution. The computational complexity of 
such algorithm is known to be 0(N 2 log TV) where N is the number of pixel 
by size in the reconstructed image. 


Theorem 2.1.2. (Central Slice Theorem) Let f be a Fourier and Radon 
transformable function. Then 

P 2 f{T cos 9, T sin 9) = T (77/) (T, 9) 

Where we mean for F 2 the 2-dimensional Fourier transform and F the 1- 
dimensional Fourier transform applied only to the variable t. 

Proof, proof can be found in [3] Ch 6, p.47; a more general result can be 
found in [3] Ch 2.1, p.10. □ 


77 

/ -- 77 / 



Figure 2.1: Scheme for Central Slice Theorem: the 2-dimensional Fourier transform 
of a function in polar coordinates is equal to the 1-dimensional Fourier transform 
made on the linear parameter of its Radon transform. 
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2.1.2 Filtered Backprojection 

Theorem 2.1.3. (Inverse Radon transform) The original function can 
he derived analytically using the filtered backprojection formula 

f= 1 -R*[R-\W\T(Rf))\ (2.1) 

where we mean that the direct and inverse Fourier transform is applied 
only to the variable t. 

Proof, see [3] Ch 6, p.49 or [4] Ch 7, p.46 for the proof. □ 

In what follows we will refer this formula as FBP. 

Given the inaccurate numerical results of this formula, we use 

w{v) = v{v)\v\ 

instead of \v\, with v a low-pass filter 1 . The reconstruction formula be¬ 
comes 

f^\n*[j^\w(y)Hnf))] ( 2 . 2 ) 

Sometimes we know w(t) = F~ l {w{v)) in order to apply the convolution 
directly to the 7 Zf, without going through the Fourier transforms, i.e. 

/ ^ ^77* [w(t) * 77/] 

An example of filter is the function w{y) = X[_l l] M with L fixed. 


2.1.3 Error bound 


The following theorem gives an error estimate for this last formula. 


Theorem 2.1.4. Let f £ C , Q O (il(0,1)), assume that the quadrature ride 
used before is exact in H' 2m , and let 77/ be reliably sampled, i.e., for some 
#€( 0 , 1 ) 

b < dm b < ir/h 


then 


f fbp = 11/ * / + ei + e2 


and 


ei| < 


1 

47T 


e*(f,b) 


i 


for a list of the low-pass filter used see the chapter about regularization. 
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M < *7(0,"OII/IU°° 

where 

£*(f,b) = \S 1 \ sup [ \J-g(a, 9)\ da 

deS 1 J\cr\>b 

and r) is any function such that 

0 < 7 /($, m) < C($) exp(—A(i?)m) 
and b > B($), B , C, A positive. 

Proof, see [3] Ch 5.1, p.104 for the proof. □ 

2.2 Iterative methods 

Given a basis of functions {bi(x)}i= i... n that interpolates the function /, i.e. 
such that 

n 

f(x) = ^2 cA(x) V.T 6 X 

t=0 

where X is properly chosen, then for the linearity of the Radon transform 

n 

nf(y) = c i n Hy) \/yeY = {(tjAi)}- 

i=0 

This is equivalent to the solution of 

Af = c 

where A(i, j) = 7 Zbi(tj, Oj) is a matrix N 2 x Ip, f is the unknown vector such 
that fi = f{xi) and c is the vector of the projection data Cj = 7 Zf(tj, 6j). 
The methods using this approach are known as Algebraic Reconstruction 
Techniques or ART. 

Using the ART approach we may have to face some problems: 

• underdetermined —> least squares 

• ill-conditioned —> regularization 

• huge —> sparseness 
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2.2.1 The Kaczmarz algorithm 

Let us choose iteratively a row k of the matrix A (usually i = k mod Ip or 
we choose randomly). Then we project the image on the vector A^. 


^(fc+l) _ x (k) 


Ci — Aj ■ x^ 

4FT 


A, 


In practice using this algorithm we reach a satisfying reconstruction only 
after having projected the image on all the rows of the matrix. In standard 
dimensions, the solution is accurate only after O(10 4 ) iteration. That is why, 
in order to increase the speed convergence, sometimes we adopt a relaxation 
parameter A, chosen properly or randomly (or with probability proportional 
to ||Tj|| 2 ) and the iterative scheme becomes 


„(*+■> =s(*>+A A-f'fA 


in this case the method is called Kaczmarz algorithm. It is shown in [4], 
[3] Ch 5.3 pl30-132 and [?] Ch 2 p.2 that the convergence is guaranteed for 
Afc G (0.2). The ART algorithm is the first method historically used, but it 
was abandoned for the FBP because of its slowness of convergence. 


2.2.2 Maximum Likelihood Expectation Maximization 

The main idea of this algorithm 2 , based on probabilistic arguments, is to 
find the x that maximizes the likelihood of the probability L(x) = P(c\x). 
This means that we are looking for the solution x which maximizes the 
conditioned probability of having our data c with an image x. 

Assuming that the noise is Poissonian on c, then the likelihood is 

lp (r*Yi 

L(x) = P(c|a;) = e -c i 

*=l Ci ' 

where c* = Ax is the exact sinogram, i.e. the projection of the exact solution 
x. Equivalently, it maximizes 


ip 

l(x) = log(L(x)) = ^2 ~{Ax)i + Ci log (Ax)i + K 
2—1 

2 see [4] for a more detailed explanation. 
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Figure 2.2: The main idea behind Kaczmarz algorithm is the projection on the 
columns of the matrix A; in the figure we used three columns as example. 


with K a constant and constraint non-negativity on x t . Kuhn-Tucker con¬ 
ditions are then: 


dl{x) 
Xj dxj 


0 if Xj > 0 


dl(x) 

dxj 


< 0 


if Xj = 0 


so that we get from the first equation 


dl(x) 


0 = Xi ~a^ = Xi 


lp Ip 

X a *'d + X 

i '=1 i =1 


N 2 

XX 

f =1 / 


';j' X j' 


or equivalently 


Xj 


x i ~ i P 


ip 

E 

i =1 


(li.jCi 


N 2 


X 1 1 X a bi 

i '=1 j '=1 
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from which comes out the iterative scheme 

( k) Ip 

(fc+i) = x j aijCj 

i i P N 2 

E“<C =1 

i'=i j'=i 

We now rewrite this scheme in a vectorial and multi-step form 

1. cf := Ax 

2. c q := c./c* (punctual division) 

3. x b = A T c q 

4- s j = Yli a i,j 

5. x'( fc+1 - ) = * x 6 ./s (product and division are made elementwise) 

Faster results are obtained by using the OSEM algorithm (Ordered Subset 
Expectation Maximization); we simply apply EM reconstruction of ordered 
subset of all the projections = 0 = (J, ©j- 

The OSEM algorithm is currently the most used in nuclear medicine 
for its speed of convergence and accuracy. 

2.2.3 Conjugate gradient 

The conjugate gradient method is widely used in all Helds for solving linear 
systems Ax = b with A square matrix of order n, symmetric and positive 
definite. In our framework the method is known as LSCG (Least Squares 
Conjugate Gradient) because we apply the conjugate gradient method to 
the normal equation of the system (??). 

In fact the matrix A T A is square, symmetric and positive definite. The 
method proposed below solves iteratively the normal equations 

A T Ax = A t c 

without having to compute explicitly A T A. The algorithm starts with an 
initialization phase: 

• given x©- 1 initial value 

• = c — Ax^ 
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Iterative reconstruction after k=1 iterations Iterative reconstruction after k=3 iterations Iterative reconstruction after k=5 iterations 



Figure 2.3: Several iterations in MLEM reconstruction of a phantom. 
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• r (o) = p ( o) _ s (o) 

• = Ap 

then begins to calculate, at each step: 

| r ( fc )|2 

• a - j^TfcTp 

• = x^ + ap^ 

• s ( fc +!) = s ( k ) _ a q( k ) 

• = A T s^ 

_ | r ( W )|2 
V | r ( fc )|2 

^ p(fc+l) — ^.(^+1) _|_ Pp( k ) 

• g(^ + l) = 

This algorithm, despite the appearance, is in most cases as fast as an 
iteration of EM. 


2.3 Kernel methods 

2.3.1 Introducing a new basis 

Now instead of the natural pixel basis we can use another basis: 

bi(x ) = B(x - Xi) = K(x,Xi ) 


with B a (essentially) compact supported and radial. 
We used the following functions: 


Function name 

/ 

Ball 

Xs(i,0)W 

Gaussian 


Wendland <£> 2,0 

(! -er)+ 

Wu ip 1.1 

(1 — £r)+(er + 2) 


with r = ||x|| 2 . 
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Figure 2.4: Original phantom and kernel reconstruction with Gaussian kernel and 
shape parameter e = 1, after 50 iterations of LSCG. 

2.3.2 Some Theorems 

Lemma 2.3.1. If fi(x) = <£>(||x||) is a radial function, then its Radon trans¬ 
form. 7 Zf is readial, i.e. it depends only on t and it is even. 

Theorem 2.3.2. If 4>(x — y) = K{x,y) is a radial function, <f> G L 1 (R d ), 
continuous, bounded and positive definite on M 2 , then its Radon transform 
7 Zf(t) is bounded and positive definite on M 1 , provided 7 Zf G L 1 (M). 

Theorem 2.3.3. (Interpolation error bound for Kernel method) 

Let f G 1)) a b-band-limited function, and let g = 7 Zf be reliably 

sampled. Let K be the interpolating Kernel function such that IZK is a 
symmetric and strictly positive definite kernel and its domain 0 be such 
that dft has regularity at least C 1 . 

Then there exist positive constants h^ and C such that, if hx.n < ho, then 

71 I AT 

11/(0 - 5 >*i( 0 IU~ < 2|5' 1 | bj-ch x , a \\nfy nK{{l) 

i =0 

with hx,n the meshsize. 
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Chapter 3 


The mathematics behind 
SPECT/CT reconstruction 


3.0.3 Radon transform and related 

We have already defined the Radon transform 




f(t cos 9 — s sin 9, t sin 9 + s cos 9) ds 


In the case of emission tomography things are more complicated. From 
Beer’s law we obtain that 


I(xi) = /(.to) exp 1 — / A(x)dx 


rx i 


'xo 


Suppose to know the attenuation coefficient, say a(x), we want to obtain 
the radioactivity /( x) by its angular projections. 

We remind the attenuated Radon transform (briefly AtRT) 


n a f(t,9)= [ e~ Va ^ e+ ^f{x) 


where T> is the Divergent beam transform defined as follows 


Vh(x,9) 



h(x + t cos 9, y + t sin 9) dt 



h(x + t9) dt 
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3.1 Analytical methods 

We recall the AtRT formula 

K a f(t,6) = f e- v «*' e+ Vf(x) (3.1) 

now, given a(x) we want to recover the value of / from its transform g = 
K a f. 

The hrst analytical solution to this problem was given in [?] by G. Novikov in 
2000. We will see the version of this formula published in [?] by F. Natterer 
in the same year. 

First of all we need the definition the Hilbert transform. 

Definition 1. Let g a suitable function, then its Hilbert transform is the 
function 

Hg(s) = 1 [ git) dt (3.2) 

it Jr s ~ t 

where the integral is meant as a Cauchy principal value. 

Note that this transform can be thought as the convolution g * h where 
h(s) = but the convolution integral does not converge. 

Now let us define the function 

h = -(/ + iT~L)lZa (3-3) 

where I is the identity operator, i the imaginary unit, and both the operator 
are meant to be applied on the linear variable of IZa. 

Now we can enunciate following Theorem 

Theorem 3.1.1. (Novikov-Natterer formula) Let g = lZ a f for f trans¬ 
formable function, and h the function (3.3). Assume a(x,y ) known, then f 
is uniquely determined by the following formula 

f(x) = ^D\ediv J e e Va i^+l\ e - h He h g) { ^ Sd) d6 (3.4) 
where S 1 = [0, 27r]. 

Proof, see [?] for the complete proof. □ 

Now let a = 0 everywhere, then g = 7 Zf, h = 0 and Va = 0 so that 

m = ^ j s Wh.0,e) de = K*W^- s g) 

which is an equivalent form 1 of (2.1). 

The proof is in [3] Ch 6 p.51. 
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3.1.1 Discretization of the AtRT inversion formula 

The discretization 2 of the (3.4) is similar to the discretization of the FPB 
formula (2.2). In fact all we have to do is to take the “attenuated back- 
projection” of the convolution product between g and the inverse Fourier 
transform of a low-pass filter. 

First of all we have to compute the function g a 

g a = 91c {e~ h/ He h g) 

Putting h = hi + ih, 2 , h\ = \lZa and I 12 = \l-LTla then 

g a = e~ hl (cos h\T-Le hl cos li 2 + sin h^Tie 111 sin ^ 2)5 

where, like in the non-attenuated case, T~i is approximated by a convolution. 
Let (Tf, a low pass filter 3 with b its cut-off frequency and E 5 = then 

U b f = * / 

once we have calculated g a we only have to compute 

f(x) = ±div [ Oe Va ^ e+ ^g a {x ■ 0, 6) d6 = K* a (g a ), 

4vr J s 1 

where the integral is computed in the sense of a discrete mean and the 
divergence operator is approximated by finite differences as explained in [?]. 


3.2 Iterative methods 

We propose a new method of reconstruction for the attenuated case. The 
main idea is the use of iterative methods as in the transmission tomography. 
We consider the outgoing Gamma-rays at angle 9 coming from a pixel with 
unit concentration of tracer. 

As we said in previous chapters, the AtRT is not equivalent to the RT of a 
weighted function. Nevertheless from the matrix A used in the transmission 
case, we can obtain a matrix B such that Be = d where B t ; j = lZ a bi(t,j, 9j), 
Ci = f{xi) the unknown term and dj = lZ a f(tj, 9j) the data. 


2 see [?] Ch 3, p.6 for more details about the discretization. 

3 again, for a list of the low-pass filter used see the chapter about regularization. 
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Attenuation phantom Activity phantom 



Figure 3.1: Analytical reconstruction of a SPECT/CT phantom data. 


Let us consider, as before, the basis b t (x) = \ P . and assume 

n 

f(x) = ^2 Cibi(x). 
i= 0 


Then for linearity 

n 

K a f(y ) = J2 dKabi{y) Vy eY = 

i=0 

Now to compute the elements of B , we have to remind that the data of 
attenuation come from a CT reconstruction, i.e., as suggested in [2], are on 
the form 

N 2 

a (x) = ]>2^Xp k (*)• 

k=1 
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For linearity its fan beam transform will be 

N 2 f+00 _ N 2 

Va{x,6) = S^g k / Xp. {x + t9) dt = Y'#*, meas{P k Pi if g ) 
k=i J[l k =i 


where = {x + tO \ s € M + } and me as is the Lebesgue measure. Now let 
us consider the AtRT of the basis 6 * 


72-A 



N 2 

E 

fc=i 


meets (P fc Fl £1" 


x,e+l 


Xp, 


This last formula is hard to compute, but it offers a useful interpretation. 
In fact if the tracer is concentrated only on the i-th pixel P u the outgoing 
rays will be subjected to an attenuation equal to the attenuated sum of the 
path length; if I in = 1, according to Beer-Lambert’s law 

/ N 2 

lout = Iin exp - ^2 9k rn,eas{P k n £f 0 ) 

\ k =l 



Now, if we consider the matrix A used in CT tomography, we can compute 
the outgoing rays from the pixel P{ in ( tj , 0 3 ) as 


Bi,j = A, exp - ^ 9k meas(P k n I tj , 0 j) 

A.j exp - ^ 9k A kl ,k 2 

where K,,) = {(AqA)} C {1,...,N 2 } 2 is the set of the indexes of the 
pixels covered by the line if. () and the index k = Ip — ((k\ — 1 )p + fo) + 1 - 
This method requires a computation time a little longer than that of the 
analytical approach because the system matrix B is calculated a posteriori. 
Moreover this matrix has often a large condition number, so the analytical 
solution seems to be preferable this time. 

Nevertheless, a regularization of the linear system can bring to a more accu¬ 
rate solution than the analytically reconstructed image. First of all we can 
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Attenuation phantom 


Activity phantom 




Reconstruction of the attenuation map 


Reconstruction of the activity 




Figure 3.2: Iterative reconstruction of a SPECT/CT phantom data with A = 0.1. 


introduce a relaxation parameter A G [0,1] in order to consider relatively 
the effect of the attenuation 

B-j = Aij exp -A ^2 9k A kl ,k 2 

Note that = A and B W = B. With this little foresight the linear 
system can lead to a more accurate solution. 

Regularization methods will be introduced in the next chapter. 
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Chapter 4 

An alternative Radon 
transform 

4.1 The Radon Transform 

In order to model correctly the tomography machine we will make some 
assumptions: 

• all rays have no width and we can consider them as straight lines. 

• all rays at angle 9 are parallel. 

• the machine runs l scans at the angles (9 \,..., 9i); the detectors move 
on a circle (typically l = 120,60). 

• for each scan j £ 1 ,... ,1 the machine gets p linear samples (typically 
p = N, or p = 128 or p = 64). 

• the total time of acquisition is relatively short so is indifferent to run 
an acquisition before or after another one. 

• the reconstructed image has dimension N x N (typically N = 128 or 
N = 64). 

Under these hypothesis and given the Beer law, we can write the Radon 
Transform 

Kf(t,6) = [ f{t9 + s6 L ) ds 

Jr 

with 9 = (cos 9, sin 9) , 9 1 - = ( — sin 0, cos 0). Now we note that this model 
is also based on another assumption: the ray is a whole straight line, so the 
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collimator “sees” all the rays in his line, also the rays behind it. 

This is not the only unreal hypothesis of the Radon transform: in fact not 
only the parallel rays will enter, because no collimator can have a resolution 
equal to zero! 


4.2 Resolution and consequences 


We write the resolution of the collimator 


Rc = D 



X + c\ 
L eff ) 


The resolution of the machine depends on the collimator resolution and on 
the intrinsic resolution (the resolution of the crystal and the electronics). 


R s =^Rl + Rj 

So we can write the resolution of the system as 



Figure 4.1: Geometrical layout of a point source acquisition by a gamma camera 
equipped with parallel hole collimator. 
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R s = \J ax 2 + bx + c 

with a,b,c some parameters we can estimate in some way. 

Now from the theory of images we see that the resolution is defined as the 
FWHM of the PSRF (Point Source Response Function), which is the data 
collected by the machine using a point source as input. 

Every image the machine sees is 

I = I 0 * PSRF 

Where Iq is the exact map and I is the map as seen by the machine. By V 
of course we mean the 2-D convolution product. 

4.3 Resolution-based Radon Transform 

For these last idea we can write the Resolution-based Radon transform as 
9) = j PRSF(x) * f{t9 + s 6 l ) ds = 

J [—r,r] 

= [[ PRSF{T,s)f((t-r)9 + s9 ± ). 

where r the the ray of acquisition. 

Now we make the (realistic) hypothesis that the PRSF function is a Gaussian 

PRSF{x) = exp(— £ 2 \x\ 2 ) 

Since we know that the FWHM of such function is FWHM = = /?(^) 

where d is the distance between source and collimator and R = R s . Hence 

rflSF(x)=exp(--th||x| 2 ) 

if we put this expression into the new transform we gets 

W(MI= /L. exp ("f?F5) m-T)s+sS±). 

since r is the distance between the parallel ray and the other parallel rays we 
are considering and the distance source-collimator is r — s. The expression 
for R 2 is simply R(d) = ad 2 + bd + c. 
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4.4 Implementation and experiments 


Since the integral is not easy to calculate for / = \ P ., we calculate the 
matrix using the following approximation 


Ri,j — 


Y ex p 

i'e[j-A,j+A] 


4 In 2 (d(j -j') 2 ] 


A 


h3 


where A is the (classical) natural pixel matrix, d is the distance between a 
ray and the next (d = N „ ), A is an integer such that 


exp 


41n2(A 2 


< 1 


and D is a matrix that contains the distances between the f — th pixel and 
the collimator when the angle is 6j. 


4.5 Experiments 
4.5.1 In silico 


Sinogram 



20 40 60 80 100 120 


Classical reconstruction 


Res-based reconstruction 


20 40 60 80 100 120 20 40 60 80 100 120 


The reconstruction using the matrix R is more accurate because the 
transformation has been done using R. 
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4.5.2 In vitro 


altezza=73 $u 128. 



altezza=73 su 128. 



In this figure we can observe that the reconstruction using the matrix R 
is more accurate, since the partial volume effect seems to be corrected. 

4.6 Conclusions 

Pros: 

• The new method is as fast as the classical one 

• From the first experiments it seems to do what we asked 

• The projection of a phantom looks like a sinogram. 

Cons: 

• New matrix as sparse as the classical one 

• There is not a strong theory that supports our choices 

• This kind of approach is completely work-in-progress. 


37 





D. Poggiali 


Numerical Methods for Medical imaging 


38 



Bibliography 


[1] Feeman, The mathematics of medical imaging: A beginners guide, Springer, 
2010 

[2] Hansen P.C., Hansen M., AIR Tools - A MATLAB Package of Algebraic 
Iterative Reconstruction Methods, Journal of Computational and Applied 
Mathematics, 2011 

[3] Natterer F., The mathematics of computerized tomography, SIAM: Society 
for Industrial and Applied Mathematic, 2001 

[4] Toft P., The Radon transform; theory and implementation, Ph.D. thesis. 
Department of Mathematical Modeling, Technical University of Denmark, 
1996 

[5] Wernick M.N., Aarsvold J.N., Emission Thomography. The fundamentals 
of PET and SPECT, Elsevier Academic press, 2004 


39 



D. Poggiali 


Numerical Methods for Medical imaging 


40 



Chapter 5 

Positron Emission 
Tomography: an 
introduction 

5.1 General introduction 

5.1.1 What is PET 

Positron emission tomography (PET) is a nuclear medical imaging tech¬ 
nique that produces an image of functional processes in the body. The PET 
machine detects pairs of gamma rays, or positrons emitted by the annihila¬ 
tion of a tracer (as the 18 FdG), introduced into the body. 



The particularity of PET stands in the possibility of recover a sequence 
of images as a time series. This allows the researchers to study the evolution 
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of the tracer kinetics, which will be treated in the next seminar. 

Since the positrons are emitted two by two, before considering a couple of 
data as the product of the same annihilation process, the machine filters 
only the pairs of positrons incoming in the same (small) temporal window 
At. This allows to reduce the effect of scattering, so the data collected after 
this filtering is quiet reliable. 

5.1.2 Mathematical modeling 

Let us introduce the mathematical models for the PET scan. 

We recall the Radon transform used in CT 

Kf{t,6)= [ f(x)dx 

Jtt,e 

and the attenuated Radon transform used in SPECT 

fcaf(t,0)= [ exp (— f a(y)dy \ f(x)dx 

Jit,e \ Jl % e (*) / 

in the PET case we receive two signals in the opposite direction i + and 
f - , belonging to the same line l. so in that case the integral of a does not 
depend on x. and then we can write 

Kaf(t, 9) = exp (- Ua(t , 9)) Kf(t , 9) 

this fact allows to correct the attenuation with ease. 

Since the image reconstruction techniques are similar to the ones used in 
CT and SPECT, in this seminar we shall only treat some post-processing 
methods as the attenuation correction and the segmentation. 

5.2 Attenuation Correction methods 

5.2.1 Soreson algorithm 

Soreson proposed a pre-correction of the data obtained from the projections. 
It is hypothesized the uniform distribution of the tracer and the uniform 
density of the tissue along each ray; this means that the algorithm will give 
out good results only when the area to reconstruct is sufficiently uniform, 
for example the brain. 

Let Pk,e and Pk.e+n two opposite projections, T the diameter body and t 
the diameter of the active region along the straight line projection. 


42 



Chapter 5. Positron Emission Tomography: an introduction 


Then, if we call N = sJPkfi • Pk,e+ir, the geometric mean of the data the 
correct value for both angles is 


N = iV 0 e 7 ? 


( t<l ) 

\sinh(/7f)y 


where 7 is the linear attenuation coefficient (e.g. the brain, if the tracer 
used and " m T we suppose 7 = 0,035cm _1 ) and / = ^ is a factor which in 
practice varies from 0.5 to 0.75. We can fix that value, for instance to 0.6 
or better to an average estimated value of ^. 



5.2.2 Chang algorithm 

The Chang algorithm makes a post-correction on the reconstructed image. 
The result of the classic reconstruction is multiplied pixel by pixel by the 
weights obtained from an estimate of the mean attenuation that the beam 
with angle 9 had been subject to at the transition from pixel (x,y). 

Then, if P is the image obtained e.g. by FBP and Po the corrected image, 
then P 0 (x, y ) = c(x, y)P(x, y) 

The correction terms are: 

C(x,y)= 

where lg i is the distance between the pixel (x, y) and the edge of the body 
at the angle and /j is the linear attenuation coefficient, supposed to be 
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constant. 

Therefore also the Chang algorithm gives good results if the tested area is 
uniform. 


5.3 Segmentation 

5.3.1 Introduction 

The aim of segmentation algorithms is to reduce an image on the form 
/ = JT \ B with B{ C M 2 some compact sets. This type of functions is 
congenial to the reconstruction because it splits the region of interest (or 
ROI) in different areas where the attenuation coefficient is constant. We can 
assign these constants if we recognize the type of tissue or we can calculate 
them. We need a segmentation algorithm that levels the image obtained by 
the image reconstruction. 


exact phantom EM reconstruction "phantomized" reconstruction 



2-norm error 36.6461 2-norm error 0.10308 


5.3.2 A simple algorithm 

The proposed “naive” algorithm is the following: 

1. P is the NxN matrix of reconstructed image, k the number of different 
levels that you want to look in [0,1] 

2. put its elements in a column vector, P = P(:) long N 2 

3. for * = 1,..., N 2 — 2 

(a) we calculate S the sum and D the difference of the elements in 
X = P(i:i + 2) 

(b) if S or D has norm less than | then P(i : i + 2) = mean(X) were 
mean is the the average 

(c) otherwise P remains the same 

(d) end of the for loop 
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4. repeat 2. and 3. for each rotation of f of the matrix P 

5.4 Level Set methods for shape recognition 

5.4.1 How Level Set methods works 

We now introduce Level Set methods, which will be useful both for the shape 
recognition and (by consequence) for the segmentation. 

As seen in [4], the main idea is to represent a closed curve of the space T as 
the zero-level set of a function T = {a? e M 2 s.t. i j>{x) = 0 }. 

We add a “time” dimension to the function in order to control the evolution 
of the curve 

T(t) = {s E M 2 s.t. (j)(x,t) = 0} 



Image courtesy of Wikipedia. 

Now for the chain rule we get 

t) ■ x = 0 

since we want the curve to expand to the normal direction n = with constant speed 

we get the differential equations 

J^ + F | V ^|=0 

with F = n ■ x a positive constant. 
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5.4.2 A Level Set method for PET 

An adaptation of the Level Set method for post-processing PET images, proposed in [2], is 
based on the MLEM reconstruction algorithm. Let the reconstruction equation be Af = c 
with / the vector of the reconstructed image values, c the raw projection data and A the 
ART matrix, as in the past seminars. 

MLEM algorithms finds the solution that minimize the functional 

L(f) = E /< " E c * MAO; + V(f) 
where V (/) is a suitable regularization term. 


Case of a unique curve 

Consider T closed curve as the zero-level-set of the signed distance function 



d(T,x) 

—d(r, x) 


x € m£(r) 

x € ext( T) 


we suppose that the image takes constant values inside and outside T, i.e. 


f(x) = kiH{<p(x)) + k 2 (1 - H(<j>(x))) 
where H is the Heaviside function so that H(<j>(x)) = Xintcr)!®) - 


Case of n closed curves 


If we have a number of closed curve it is a bit more difficult. For i = 1 ,,n, let the 
binary representation of * — 1 be bin(i — 1) = (6j, ..., bh). Then / can be written as 


2 71 n 

f( x ) =n fc< ri- Ri (A( a; )) 

i =1 i --1 


where 


Ri((f)j) — 


H(4>j) if b) = 0 
1 -Hifa) if bj = 1 


We calculate with the chain rule Sr ~. 

We easily see that = e — A T (b./Af) plus the derivative of the regularization term. As 
regularization term we use the length of the curves V{<j>j) = f \VH(0j)\ whose derivative 
is known 


dV 

d<t>i 


= - V- 


( V& 


VI 




With these tools we can compute JV- , and use it in the following iterative algorithm: 


1. Choose initial level sets (j >\ 01 and ’time step’ At. 

2. Update the level sets 4>\ n+1 ' > = </>A — /\tAA 

3. Reinitialize the level sets. 
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Figure 5.1: An example of shape recognition using the level set approach. 
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Chapter 6 

Kinetics of the Tracer in 
PET 

6.1 What are compartment models 

A (multi-)compartment model is a type of mathematical model used for describing 
how materials or energies are transmitted among the compartments of a system. Each 
compartment is assumed to be a homogeneous entity. Hence a compartment represents, 
in medical application, a group of organs exchanging material at the same speed one with 
another. 

The kinetics of the material is described by an ordinary differential equations system 
with as unknown functions the concentration of the considered material in each compart¬ 
ment. 


6.1.1 A simple example 



Qi/v 


Figure 6.1: An example of compartment model. 


51 




D. Poggiali 


Numerical Methods for Medical imaging 


The equations for the model of Fig.6.1 is 


Q 1 = Pi — K 01 Q 1 + D8(t) 

Ci = sm 


Q i (0) = Qio 


6.1.2 General model 

The general model for compartment model is the following: 


Qi = Y, - I] F a(Q) + p i + u i(t) 

i¥=i 

Ci = or Ci = 0 


Qi( 0) — QiO 


6.2 Tracers 

6.2.1 Definition 

A tracer is a small quantity of substance used for tracking purposes. 
The tracer: 


1. can be measured independently from the traced 

2. follows the dynamics of the traced and do not perturb the state of the system (small 
quantity) 

3. is kinetically indistinguishable from the traced 


6.2.2 Indistinguishably principle 

From the third assumption we can get the following principle: 

P (a particle leaving j is tracer) = P (a particle in j is tracer) 
this means that 


h 


1j 


equivalently 


Fij + fij Qj + qj 




We can write the general model for the kinetics of the tracer as follows: 
f 4i = ~ + w *( t ) ’ <?;(0) = <7;o 


r — n(t) 
— Vi 


j^i 

a = 0 


If we suppose that all the functions kij(Q) = kij are constant then fij(Q) = kijqj and 
the model becomes: 


qi = Y kgqj — Y kjiqj + Uj(t) , qi(0) = qm 


Ci 


j^i 

_ 


Ci = 0 
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Figure 6.2: A general model for tracer kinetics. 


6.2.3 A more general view 

Under this hypothesis the model goes under the category of Linear Time-Invariant model 
in state-space form LTI state-space 

f x = Ax + Bu , *(0) = xo 

\ y = Cx 

Where u is the input, y the output and the system is determined by (A, B, C). 
Theorem 6.2.1. The system (A,B,C ) is equivalent to the direct input-output relation 
y(y) = Cx(t) = (C exp(tA) B ) * u(t) 

Proof. If we consider x, the Fourier transform of x we get: 

sx = Ax + Bu 


which is equivalent to 

x = (si — A) -1 Bu. 

and finally 

y = Cx = C (si- A)' 1 B 

that leads to the thesis. □ 

So we define the impulse response function 

4>(t) = C exp(tA) B 

and its Fourier transform, the transfer function 

4>(s) = C (si-A)- 1 B 

so that y = (j>*u and y = cpu. 
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Figure 6.3: Is this model identifiable? 


Moreover under the hypothesis the the input is a bolo 

u(t ) = DS(t) 


we get that 


y(t ) = (j>{t) * u(t) = * D5(t) = 


6.3 Identifiability and identification 

6.3.1 Some definitions 

Definition 2. The model (A,B,C) 

( x = Ax + Bu , x(0) = xo 

\ y = Cx 

is (uniquely) identifiable if it is possible to determine (uniquely) the nonnegative entries 
of A, given B,C, and the exact data xo,y,u. 


Definition 3. 1. A Compartment k is in/out-put reachable if 3 a path from an in- 

out-put compartment to k. 

2. A system ( A,B,C ) is in/out-put reachable if every compartment is in/out-put 
reachable. 

Example 1. the system of Fig. 6.4 is input reachable but not output reachable. 
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6.3.2 Necessary and sufficient conditions 

Theorem 6.3.1. (Necessary conditions): Let(A,B,C) be an identifiable system, then: 

1. The system is input reachable 

2. every compartment with an outgoing path is output reachable 


Example 2. Example: the system of Fig.6.4 could be reachable. 


Observation: The transfer function <f>(s) is a rational function of the form 


fi(s) 


with N £ V n ~ 1 ,D £ V 

D(s) 


where n = dim(A). 
In fact 


fi(s) =C {si- A)- 1 B = 


det(sl — A) 


n — 1 

and adj(A) = RkS k 
k = o 

Theorem (sufficient conditions 1) If the equation of the numerator and denominator 
have a (unique) solution, then the system (A,B,C) is (uniquely) identifiable. Theo¬ 
rem (sufficient conditions 2) If the equations (o) = CA k B , forfc = 0, ... ,n — 1 
have a (unique) solution and the denominator is never null, then the system ( A,B,C ) is 
(uniquely) identifiable. 


Example 3. The system of Fig.6.4 

/ 0 ki k 3 \ 

A= 0 ~{ki + k 2 ) 0 

V 0 k 2 0 / 

B = (0,1,0), C = (0, 0, 1 /V 3 ) 
is identifiable, but not uniquely. 


6.3.3 Identification 


Identification can be made using the results of the two sufficient conditions, but such 
techniques has revealed to be numerically inaccurate. So the calculation of the parameters 
is given by a (weighted) nonlinear least squares method. 


A = argmin ||j/ — j/a|| 2 = \\y~4>A*u\ 


= ||y — (C exp (tA) B ) * m| 


6.4 Sokoloff Model 

The Sokoloff model (1970s) is the main model studying the kinetics of the tracer 18 FdG 
in the brain. It uses a sequence of PET reconstructed images as input. 
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Figure 6.4: Scheme for the passage from PET data to time-dependent data. 


C = V b C p + (l-V b )[C s + Cjb\ 



The model is the following 

/ —fci k 2 0 \ 

A = ( ki — {k 2 + kz) 0 | , 

V 0 k 3 0 / 

B = (1,0,0), 

C = (14,1 - Vi, 1 - V b ). 


Is this model identifiable? 


det(sl — A) = s 3 + (hi + k 2 + k 3 )s 2 — k\k 2 
C adj(si — A) B = V b s 2 + {V b ki + k 2 )s. 

Answer: yes it is. 

From the values ki, k 2 , k 3 we can compute the value of the local glucose metabolism 


mr = k t? 


where with this model K = , C p is concentration of glucose in the plasma (in 

stationary state) and LP is the lumped constant. 
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Chapter 7 

An introduction to MRI for 
mathematicians 

What is MRI 

Magnetic Resonance Imaging (MRI), elderly known as Nuclear Magnetic Resonance (nMR 1 ) 
is an in vivo medical imaging system that uses the resonance of magnetic fields. This exam 
is non-invasive and much secure compared to the analog X-ray tomography. In fact the 
electromagnetic rays in use have a maximum frequency of about 10 5 Hz (which is in the 
order of radio waves), which is much less than the X-ray frequency, about 10 16 Hz. A 
higher frequency means a shorter wavelength, a greater energy and by consequence a 
greater mutagen power. When we say that RMI is “safe” we mean that the probability of 
a genetic mutation of the RNA a cell is low. To make more safe the using of X, 7 and f) 
rays we have to reduce the time of exposition. 

7.1 A Physical introduction 

7.1.1 The spin 

The nucleus of an atom has a property called spin. This means that we can consider 
an atom as a spinning wheel. This property leads also to an electromagnetic effect: each 
nucleus generates a small magnetic field. This means that we can also think of nuclei as 
a magnet with North and South poles on the rotation axes of the nuclei. 

For our purpose we will consider the nuclei of hydrogen (H), because is has an odd number 
of protons and so the magnetic moment is not null and also because is very common. 

So a proton p has an intrinsic angular momentum J p and a magnetic momentum /r p . 
There exists a direct proportion between these two quantities 

}i p — 7 /' Jp 

where "y p is a constant called gyro-magnetic ratio and it is related to the kind of atom in 
examination. For instance the hydrogen has a ratio equal to 7 h — 2.68T~ 1 s~ 1 . To be 

1 the ‘n’ has been remover to underline the complete safety of this kind of exam. 
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more precise the gyro-magnetic ratio depends on the chemical bound of the atom 


7p = (1 - 

where the constant a £ [ 10 ~ 6 , 10 “ 4 ] is called chemical shift. 

Now if we turn on a stationary and uniform magnetic field Bo, the quantum mechanic 
expected value of the intrinsic angular momentum will satisfy the law 


d < J p > 
dt 


=< n P > xB 0 


so, given the proportion we obtain 


d < fj, p > 

dt 


7 < Up > x B 0 


since we can not measure the momentum of every single atom, we define the (total) 
magnetic momentum M as 

M = 'y ' fj, p 

p 

and then we can formulate the differential equation 


M = 7 M x Bo 


(7.1) 


which will be fundamental for our purposes. Even if the nuclei are under a quantum effect, 



Figure 7.1: A nucleus of hydrogen can be seen as a magnetic wheel, 

a classical approach will be more feasible and will not lead to great errors. 

7.1.2 Magnetic fields and relaxation 

In absence of an external magnetic field the nuclei arrange themselves randomly and their 
momentums mutually annul, so M = 0. Since the magnetic field of the Earth is weak 
(from 2 to 7 -ICE 5 T), we can consider that a normal body has a magnetic momentum 

M « 0. 

Otherwise, if we turn on a stationary magnetic field along the z-axis Bo = (0, 0, Bo), the 
nuclei will show a sort of herd mentality. In fact their rotation axes arrange themselves 
parallel to the axis 2 . Some of them align themselves with the same verse of B 0 (low 
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0 0 0 C 1 

v 0 


Figure 7.2: Herd mentality of the nuclei under a magnetic field. 


energy), and some others will take the opposite direction (high energy). As we will see in 
the next section, the most part of the nuclei will arrange in a low energy position. 

A more accurate description of the motion reveals that the axes has a precession 
along the direction of the external field at the same frequency u>o, called the Larmor 
frequency, but with different phases. 

The precession tends to the equilibrium (0, 0, M eq ) exponentially, with a recovery time Ti 
(spin-lattice) along the axis z, and with a loss time T 2 (spin-spin) along the plane xy. 
These two times are different because of the different phase of the nuclei. 

The parameters Ti and T 2 depends from the chemical bound of the Hydrogen molecules, 
then from the tissue. 

The main problem we will face is that the precession is not visible for the total momentum 
because the convergence to the equilibrium is too fast. 


7.1.3 Attempt for seeing the relaxation 

The disposition of the nuclei in the verse of the magnetic field causes the precession motion 
to be difficultly measurable. To make the precession visible we have to study deeper, the 
alignment to the axis of the nuclei. According to the Boltzmann factor the probability 
that a nucleus arranges himself along the 2 -axis is equal to the negative exponential of 
the ratio between magnetic and thermal energy, 

r>{ rp \ ( Emagn \ ( T-Bo 

PiE,) = ’ xv {-EZZ ) =exp (-55iJ 
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where H = 1.05 • 10 -34 Js is the normalized Plank constant, kb — 1.38 • 10 -23 JK ~ 3 is the 
Boltzmann constant and Tk is the absolute temperature in Kelvin degrees. 

Let us calculate this value for a typical MRI experiment, Bo = IT, T = 300A'. 


P(Ei) = exp 


/ 3 s _ 1 T _1 -IT- 10 -34 Js \ 

V 1 ■ 10~ 23 JK- 1 ■ 3 ■ 10 2 A' ) 


exp (— 10 ~ 13 ) 


for the first order Taylor approximation we can compute 

P(Ei) 1 - 1CT 13 . 

This means that only one atom over 10 13 chooses a high energy position, and the magnetic 
momentum tends quickly to the 2 -axis. 

From the Boltzmann factor comes an idea to make more visible the precession: we can 
decrease the room temperature or increase the external magnetic field. 

For having a visible change we have to move two parameters in values that are physi¬ 
cally unacceptable. In conclusion we have to find another way to make the precession 
measurable. 


7.1.4 RF pulse and resonance 

To move down the magnetic momentum we have to resort to a Radio-frequency pulse (RF- 
pulse), a second magnetic field JBi of a minor intensity than the static field and rotating 
along the plane xy at exactly the Larmor frequency u>q. In this way a physical phenomenon 
will occur, the magnetic resonance. In fact a huge number of nuclei will gain energy 
and and the magnetic momentum M will separate from the 2 -axis. 

Moreover with the RF field turned on the phase of the precession will align, while when 
we turn off the RF field there will be “dephasing”. 


7.2 A Toy-model: the Bloch equations 

To see mathematically phenomenons described in the previous section, we need a toy- 
model, a set of mathematical equations that describes the motion of the aggregate mag¬ 
netic momentum M = (M x , M y , M z ). So let us consider the equation (7.1) applied to an 
external field B and study the solution of these equation for several external fields. 


7.2.1 Static field 

Let us consider the case that only the static field along the axis 2 is on. In this case we 
set B = (0, 0, .Bo). If we add a relaxation term that describes the polarization along the 
axis 2 (7\), and the decay along the au/-plane (T 2 ) we get the Bloch equation 


M = 7 M x B 


( M x My M z - M eq \ 
V T 2 ’ T2 ’ Ti ) 


(Bloch-1) 


we can write this equation component-wise 

f M x = M y B 0 -^ 

i M y = M X B 0 -f§ (7.2) 

[ M z = 
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The solution to this ODE is 


M x 

= e t / 1 ' 2 (JMj,0 cos cuot — M y , 0 sin cuo t) 


My 

= e~ t / T2 (M x ,0 sin wo t + M y , 0 coscuot) 

(7.3) 

M z 

= M Zf 0 e- t/T i + M eq ( 1 - e~ t/Tl ) 



where cuo = — 7 -Bo is the Larmor frequency. This is clearly a precession along the 2 —axis at 
the Larmor frequency, and the solution tends to the equilibrium (0, 0, M eq ) in a exponential 
way, so the precession is very hard to see experimentally. 


7.2.2 RF pulse 

We turn on a RF pulse, rotating with angular speed cu in the plane xy, we now have 
B = (B 1 coscut, B 1 sin cut, Bo). Since we consider what happens while the pulse is on, we 
can ignore the relaxation terms and write 


M = 7 M x B 

again, we rewrite the equation element-by-element 

{ M x — 7 (M y Bo — M Z B\ sincut) 

M y — 7 (M z Bi coscut — M x Bo) 

M z = 7 ( M X B\ sin cut — M Z B\ cos cut) 


(7.4) 


(7.5) 


now, for simplicity, we study the system in a reference frame rotating along the axis 2 like 
the RF pulse. 


el = 

(coscut, 

sin cut, 

0 ) 


e2 = 

(— sincut, 

coscut, 

0 ) 

(7.6) 

e3 = 

( 0 , 

0, 

1 ) 



Now if we consider the functions 

u(t) 
v(t) 


— Rujt 


M x (t ) 
My(t) 


(7.7) 


where R^t is the matrix of rotation of angle cut. 

Substituting the functions (7.7) in the system (7.5) we get 

{ u = (-/B 0 + c o)v 

v = — ( 7 B 0 + cu) u + -/B 1 M Z 

M, = —-/Biv 


(7.8) 


now we observe that if the angular speed of the RF pulse is exactly the Larmor frequency 
cu = cuo = — 7 Bo the system is simplified and we obtain the second version of Bloch 
equation 

f u =0 

< i) = 7 B±M Z (Bloch-2) 

[ M z = —-/Biv 

this means that the system is affected by resonance, since and the aggregate momentum 
tends to go move down from the equilibrium position. This result is known as Larmor 
Theorem. The solution of this equation is simple 


u = Uo 

v = Vo coscuif — M 2 ,o sincuif 
M z = vo sincuit + M z n coscuit 


(7.9) 
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Figure 7.4: The solution of the second Bloch equation in the non-inertial reference 
frame. 

where, as one can guess, c o g = —yB g . 


7.3 Pulse sequences 


As we can see from (7.9), an RF pulse has the effect to flip the vector M from the axis 2 
of an angle 9 = uit after a time r. Thus we can apply following pulses: 

• a general 9- pulse: t = 

• a f-Pulse: n = 

After this pulse the aggregate magnetic momentum is 

M(n) = {uo,M z ,o,vo). 

• a 7r-pulse: T 2 = ^ . 

After this pulse the aggregate magnetic momentum is 

M(t 2 ) = («o, -vo, 

We can play with pulses to get more informations about the precession. The following 
sequences are the most popular ones: 
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1 . inversion recovery 

2 . saturation recovery 

3. spin echo 

7.3.1 Inversion recovery 


Inversion recovery: after the we apply a 7 r-pulse, we wait a time r, then we apply a 
^-pulse and we measure M. We wait for the system to get back to the equilibrium of 
(7.3) and we repeat several times to increase the signal-to-noise ratio and for several values 
of t. This sequence is shorted as follows 

(7T — T — ^ — AT — too)n 

where AT stands for acquisition time, and t a0 is a long time needed to get back to the 
equilibrium, usually too « 4Ti. 

From the third component of (7.3) we get 

M z (r) - M eq = (M z ,o - M eq ) exp(-T^). 

Ji 

This means that we can obtain — 1/Ti as the slope of log— M eq ). 

7.3.2 Saturation recovery 

Saturation recovery: This sequence is shorted as 

{ l_ Hs _ T _l_ AT _Hs) n 

Where Hs is a pulse that destroys the homogeneity of the field Bq\ after a time r the 
magnetization is partially recovered, the pulse fplips it on the xy-plane. Again, we obtain 
— 1/Ti as the slope of log(M z (r) — M eq ). 


7.3.3 Spin echo 

Spin echo: This sequence is shorted as 

(| - t-7 r )„ 

After 2 t the amplitude of the transverse magnetization will be relaxed by a factor 

exp(2r/T 2 ) 

T2* is the real time constant of loss in the xy-plane, because the field Bo is not really 
homogeneous. The relation between T2 and T2* is quite simple: 

1 1 AT, 

- =-h 7 A Bo 

T2* T2 

where A B 0 is the maximum variance of B 0 - 
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7.4 K-space and reconstruction 

7.4.1 From Bloch equations to imaging 

The aim of reconstruction is to create a 2D map of the proton density p(x) on each slice. 
For this purpose we consider only the transverse component of M, 

M* = M x + iM y 

and same thing for the signal received from the coils 


S - S X T f'Sy 

For Faraday’s law the magnetic field generated by the body induces an electromotive 
force to the coils 


S(t) = ^loo P - k Jtl oop p{x)Mdx 

Now, if we consider a tiny slice over the z-axis, we can rewrite the previous formula as 

s(t) = ki^- f p(z)M*(z)dz 
dt 7loop’ 


with zeC . 


If we apply a Il-pulse, then a gradient for the selection of a slice, we get 


then we obtain 


M* = exp ( — — + iojot + iojgt 


s(t) = fc 2 exp ( — — _|_ ) / p(z) exp(iw g t)dz. 

\ 1 2 / J loop’ 


We remember that uj g = —jBg = —yG ■ x then, if we set 


G — G x + iG y 


we get 


s(t) 


&3 exp 


\--l—+iuJat\ / p(z) exp(— it'yG ■ z)dz 
\ 7 2 J Jloop’ 


7.4.2 K-space 

If we look at the previous formula, we have found that 


where of course k(t) = —P-G. 
So we store the values of 




exp(A 

2-D IFFT algorithm to recover the slice image. 


s(t) oc Tl [/E»](fe(t)) 

in a matrix called K-space and we apply the 
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Reconstructed image 




50 100 150 200 250 


50 100 150 200 250 


An example of K-space and reconstructed image. 


Simulated raw data 2DFFT Reconstructed image 



50 100 150 200 250 20 40 60 80 100 120 20 40 60 80 100 120 


We can also use a phantom to obtain a simulated K-space. 

7.5 TR, TE and Sequences 

7.5.1 Parameters 

MRI is a multi-parametric imaging system: in fact you can play with some parameters to 
obtain different contrast between the tissues. 

We now introduce some of these parameters. 

TR, Repetition Time is the time interval between the application of a pulse and 
another. 

After TR the longitudinal magnetization of the vector M will be 

M z {t) = M eq (l-e- TR / T1 ^ 
so each FID signal will be proportional to 

(i-e -: TR/T1 ) 

This means that the longitudinal magnetization will not reach M eq unless we set TR ^$> 
4T1 

TE stands for echo delay time (or time to echo). Instead of making the measurement 
immediately after the RF pulse (impossible), we wait a short period of time and then 


68 








Chapter 7. An introduction to MRI for mathematicians 


make the measurement. 



90° 180° 

RP___/T\_ 

90° 

_ u _iT\_ 

--E 

h- 

: e 

1 i 

1 TE/2 | 

Echo 1 

oignai j l '\JI 


Gx—- 1 1 - 


A complete example of spin-echo sequence. 


7.5.2 Sequences 

• Tl-weighted short TR 

Useful in scanning brain to see contrast between gray and white matter. 


• T2-weighted long TE, long TR 

Used to see the contrast between water and fat. 


• T2*-weighted long TE, long TR 

We use it to see particulars in venous blood. 

• Proton density weighted short TE, long TR 

We try not to consider the effect of T 1 or T 2. 
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TR 

Time -> 


h 2 0 

Fat 

Solid 
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Chapter 8 

Dual-Modality Imaging 


8.1 Introduction 


PET/MRI is a new, arising hybrid-imaging technology that incorporates MRI morpho¬ 
logical imaging and PET functional imaging. 

At the moment 3 companies declare availability of their PET/MRI imaging devices: 


• PHILIPS 

• SIEMENS 

• GE 
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8.2 Reconstruction in PET/RM 
8.2.1 Classic MLEM 

Consider the linear system Ax = b, with b projection (PET) data and x the image to 
reconstruct Maximum Likelihood Expectation Maximization (MLEM) is based on 
a probabilistic argument (the noise is assumed to be Poissonian). 

lP 

L(x) = P(b\x) = -h~ e ~ bi 

Oi ! 

t=l 

where b* = Ax is the exact sinogram, i.e. the projection of the exact solution x. Equiva¬ 
lently, it maximizes 


ip 

l(x) = log (L(x)) = ~(Ax)i + bi \og(Ax)i 

i= 1 


Kuhn-Tucker conditions are then: 

dl(x) 


= 0 if £, > 0 


^<0 if Xj = 0 

dxj 

From the first equation (the second is always satisfied) we get 

/ \ 


dl(x) 

0 = Xj — = Xj 

dxj 


this leads to the iterative scheme: 


ip 


ip 


- J2 °i'.i + I] ■ 

i'=1 i =1 


Q*i, j bi 


y! a iyj 


j '=1 


Ip 




bi 


E -'j ’ E <•« 


j '=i 


If we write the scheme in a vectorial and multi-step form 

1. b f := Ax {k) 

2. b q := b./b* (punctual division) 

3. x b = A T b q 
4- s i = Ei a «.j 

5. a4 fc+1) = x ( k K * x b ./s (product and division are made elementwise) 
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8.2.2 Minimum Cross Entropy 

The goal of a reconstruction iterative scheme is to minimize the distance D(b,Ax). For 
instance the least-squares method proposes to minimize the euclidean distance 

D] (b, Ax) = Y ~ b i 

i \ .7 

Note that using this distance leads to the normal equations. 

The Bayesian methods are similar, but they add an extra penalty term D 2 to the 
objective function 

J 3 (x) = Di ( b, Ax) + /3D 2 ( x) 
with the constant P > 0 set up by the user. 

Sometimes D 2 can be the distance between the current image and a prior image model 
p. In this case the objective function is 

Jp(x) = D x (b,Ax) + /3D 2 (x,p). 

The cross-entropy or Kullback-Leiber distance is defined as it follows 
S(U, V) = Y Ui l°g (Ui/Vi) — Ui+ Vi. 



Note that minimizing S(b, Ax) is equivalent to minimizing the log-likelihood function l[x) 
described above. 

So we minimize 

J 3 (x) = S ( b , Ax) + f3S(x,p). 

Minimizing this function leads to the iterative scheme: 


Y a i' 


Y - p log {Xj/Pj 

i / , a i,j' x j' 


This scheme is called MXE1 and converges provided the non-negativity Kuhn-Tucker 
conditions 

bi 


E 


Y! 


> P log (Xj/pj 


so for p too large this condition is violated. 

Another iterative scheme for solving the optimization problem is 


Xj Pj exp 


-Me-o-e 


O'i ,j bj 


'Yhj' 


this scheme is called MXE2. MXE2 satisfies the non-negativity constraint for all p > 0. 
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